Gravitational Waves From the End of Inflation: Computational Strategies 



Richard EastheiQ and John T. Gibhn, JiQ 
Dept. of Physics, Yale University, New Haven CT 06520 

Eugene A. Linf] 

Department of Physics and ISCAP, Columbia University, New York NY 10027 

(Dated: February 4, 2008) 

Parametric resonance or preheating is a plausible mechanism for bringing about the transition 
between the inflationary phase and a hot, radiation dominated universe. This epoch results in the 
rapid production of heavy particles far from thermal equilibrium and could source a significant 
stochastic background of gravitational radiation. Here, we present a numerical algorithm for com- 
puting the contemporary power spectrum of gravity waves generated in this post-inflationary phase 
transition for a large class of scalar-field driven inflationary models. We explicitly calculate this 
spectrum for both quartic and quadratic models of chaotic inflation, and low-scale hybrid models. 
In particular, we consider hybrid models with an "inverted" potential. These models have a very 
short and intense period of resonance which is qualitatively different from previous examples studied 
in this context, but we flnd that they lead to a similar spectrum of gravitational radiation. 



I. INTRODUCTION 

The universe has been transparent to gravitational ra- 
diation since its energy density dropped below the Planck 
scale. Since they are unhindered by optical opacity, grav- 
itational waves from the early universe can constitute 
a stochastic Gravitational Wave Background, or GWB, 
that complements the familiar Cosmic Microwave Back- 
ground (CMB). Like any other form of radiation, the 
number density of gravitons drops as the cube of the 
scale factor, a, while their energy is redshifted by a fur- 
ther factor of a, leading to pgw ~ 1/a*. If a remnant 
background of gravitational radiation is laid down near 
the Planck scale, this would occur long before the end 
(or even the onset) of inflation. Consequently, it would 
now be radically diluted, and thus entirely unobservable. 

On the other side of the ledger, a number of processes 
in the early universe can generate gravitational waves, 
yielding a GWB which is potentially detectable in the 
present epoch. Unlike the CMB, the GWB is strongly 
model dependent, and there is no consensus as to its 
likely form. Since the existence of the CMB is a generic 
prediction of hot big bang cosmology, the detection of 
black body microwave radiation by Penzias and Wilson 
[1 played a key role in establishing the big bang as the 
dominant cosmological paradigm. By contrast, the ab- 
sence of clear expectations for the GWB implies that a 
stochastic gravitational wave background which is not as- 
sociated with "standard" astrophysical sources will pro- 
vide a sensitive tool for discriminating between different 
models of the early universe. 

The existence of gravitational radiation is inferred from 
studies of binary pulsar systems, whose orbits decay at 



a rate entirely consistent with the energy they are ex- 
pected to emit as gravitational radiation [2 . A similar 
inference would be drawn if a B-mode polarization contri- 
bution was seen in the foreground-subtracted CMB sky. 
However, in neither case is one directly observing the 
distortions of space caused by gravitational waves propa- 
gating through the universe. This is the province of "di- 
rect detection" experiments such as the currently oper- 
ating Laser Interferometer Gravitational Wave Observa- 
tory (LIGO), future space-based observatories, like Laser 
Interferometer Space Antenna (LISA) and Big Bang Ob- 
server (BBO), or detectors which are sensitive to the de- 
formations of solid objects or cavities. The frequency 
range of these experiments is largely determined by their 
physical size.^ Finally, proposals such as the Parkes Pul- 
sar Timing Array represent a hybrid appoach, using an 
array of pulsars as a "mesh" of sensitive clocks which 
effectively measure changes in their distance induced by 
a gravitational wave background. This "instrument" is 
sensitive to gravitational waves in the nano-Hertz range, 
or wavelenghts of O(IOO) hght years [5]. 

The best known candidate source for the GWB is the 
quantum mechanical fluctuations of spacetime during in- 
flation, which provide a near scale-invariant spectrum of 
gravitational waves in the present epoch. The existence 
and near scale-invariance of this spectrum is a key predic- 
tion of inflation. However, the amplitude of this spectrum 
is proportional to the inflationary energy scale, which is 
a very poorly constrained parameter, as there is no con- 
sensus or "standard model" of inflation. For convenience, 
this amplitude is expressed via the tensor:scalar ratio, or 
r - the denominator is known to be 10~^ from mea- 
surements of the CMB and large scale structure. Current 



1richard.easther@yale.edu| 
liohn.siblin@yale.ed u 
Teugene . a. limijigmail . com 



^ LIGO is an exception here, as it operates in a Fabry-Perot mode, 
so its characteristic wavelength is actually some 10^ times larger 
than the physical size of its interferometers. 
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data for the CMB B-mode suggests that r < .3, which 
already rules out well-known models of inflation [H [5]. 
There are two contending theoretical arguments as to 
the most likely value of r. Insisting on an inflationary 
potential which has a stable minimum with a vanishing 
potential energy (so as not to contribute a cosmologi- 
cal constant) and a simple algebraic form suggests that 
r 0.001 [51. In practice, this limit depends very much 
on one's definition of simplicity.^ Conversely, attempts 
to build inflationary models inside string theory typically 
yield r < 10~^°, which can be traced to the need to avoid 
a trans-Planckian vev for the inflaton [71 |H1 IS] • In this 
case the inflationary GWB is unobservable by any pro- 
posed experiment. Consequently, the failure of a BBO 
style experiment to detect a near scale-invariant back- 
ground would rule out many popular models of inflation, 
but would not undermine the overall paradigm. 

In this paper, we consider the generation of gravita- 
tional waves at the end of inflation [TUl dU d^]. Iron- 
ically, in models with a low inflationary scale, it turns 
out that this signal - if it exists ~ can naturally fall into 
the frequency range probed by BBO, and other future 
experiments including LISA and an upgraded LIGO. By 
contrast, the conventional inflationary spectrum is effec- 
tively invisible in these scenarios. Consequently, this sig- 
nal extends the range of models which experiments like 
BBO would be able to test. Moreover, since this sig- 
nal depends on the details of the mechanism by which 
inflation ends, it provides a window into an otherwise 
unobservable epoch in the early universe. 

In most models of inflation, the energy of the uni- 
verse is dominated by the potential energy of a homoge- 
nous and isotropic scalar field, (j). As inflation ends, this 
energy must somehow be used to generate elementary 
particles and reheat the universe. The specific mecha- 
nism by which reheating occurs is strongly model depen- 
dent. Originally, the creation of other particle species was 
thought to occur slowly, since the inflaton field is neces- 
sarily weakly coupled to other fields. Although the bot- 
tom of this potential has a vanishing energy, the field(s) 
oscillate homogeneously, with their motion damped only 
by the expansion of the universe. Parametric resonance 
provides an efficient mechanism for extracting energy 
from these oscillations, and converting it into excitations 



^ The analysis of 6 concentrates upon potentials which are no 
more than quartie in the inflaton Models with a low scale 
are typically of the form ~ Vb — Ai/)* with any quadratic term 
strongly suppressed [B]- In this case, one needs a i})^ term in order 
to construct a stable minimum. Since V{(t>) is the effective poten- 
tial it can contain higher order terms and/or logarithms. This al- 
lows one to explicitly construct low-scale single field models with- 
out including non-renormalizable terms in the action, although 
one might object that loop contributions will only be dominant 
in models where the tree-level terms have been carefully tuned. 
However, high scale models of inflation will be severely tested by 
the next generation of CMB polarization experiments, and this 
question will ultimately be decided by experiment. 



of other fundamental fields [T5^. This occurs even when 
the tree-level couplings between the inflaton and other 
fields are very weak, and the amplitudes of the momen- 
tum modes are approximately described by Mathieu-like 
functions. 

In the classical picture, the momentum modes are sim- 
ply the Fourier components of the corresponding fleld - 
and if these amplitudes are growing exponentially, it fol- 
lows that the fields and their associated energy densities 
will become increasingly inhomogeneous during the reso- 
nance phase. This leads to the emission of gravitational 
radiation, and this process forms the topic of our inves- 
tigations here. The first connection between these two 
issues was made by Khlebnikov and Tkachev [10], who 
predicted a gravitational wave signal whose peak ampli- 
tude was approximately Q,gw ~ 10^^'^ for quartie, A^"* in- 
fiationaxy models, at (present day) frequencies of around 
1 GHz. This amplitude is substantial, but as the required 
strain sensitivity needed to detect a signal of fixed am- 
plitude scales as the cube of the frequency, its detection 
presents an enormous technological challenge. 

In early 2006, Easther and Lim revived this topic, and 
argued that the amplitude of any preheating signal could 
be essentially independent of the inflationary scale, while 
its (present day) frequency was proportional to the en- 
ergy scale of inflation [T^]. Thus, while a GUT scale 
model would be peaked near GHz scales, resonance fol- 
lowing inflation at 10^ GeV would lead to a signal near 
the LIGO band, while lowers scales would overlap with 
the proposed range of BBO. In [T^, this proposal was in- 
vestigated in a limited way by considering models where 
the inflationary scale differed by a factor of \/T0 using 
a fairly rudimentary numerical algorithm, derived from 
that employed in |10j . This algorithm had a number of 
limitations, in that it did not allow the gravitational wave 
signal to be computed at arbitrary times, as it considers 
the power generated in a series of four dimensional space- 
time "boxes" , and was based on a formula for the radi- 
ated power in gravitational waves which is only strictly 
valid in a non-expanding universe. In the discussion be- 
low, we denote this the "box" algorithm. 

In [13], the present authors introduced a new algorithm 
which directly evolved the tensor component of the met- 
ric perturbation, and explicitly conflrmed that a preheat- 
ing signal could be visible with future versions of LIGO or 
BBO when the inflationary scale is low enough, and that 
its amplitude was essentially independent of the inflation- 
ary scale. This algorithm directly solves for the evolution 
of the momentum-space metric perturbations sourced by 
the transverse-traceless part of the stress-energy tensor, 
and we refer to it as the "spectral method." We assume 
a Friedman-Robertson- Walker background metric with 
perturbations in synchronous gauge, 

ds^ = dt^ - a^it) [5,j + hij\ dx'dx^ . (1) 

where h^j is Transverse-Traceless or, 

hi = 0, hi, = 0. (2) 



3 



The analysis here assumes that the universe has only 
scalar field components, but generalizes to arbitrary T^i^ 
- including other situations with substantial inhomogene- 
ity, such as phase transitions, bubble collisions or turbu- 
lence. In numerical analysis, "spectral methods" refer to 
algorithms that decompose the quantity of interest into 
a set of basis functions, and solve for their coefficients to 
describe the system as it evolves. In our analysis, we ac- 
tually evolve the fields in position space in an expanding 
background, using LatticeEasy [15]. We then use this 
solution to source evolution of the metric perturbations, 
so the overall approach is essentially a hybrid scheme. 

Two further numerical algorithms have been developed 
to compute this signal - a welcome development, given 
the intricacy of the relevant calculations. In [TB] , Garcia- 
Bellido and Figueroa restated the scaling argument of 
[12j and computed the gravitational waves sourced by 
bubble-collisions at a variety of scales. The amplitude 
of the signal seen in [T^ is not independent of the in- 
flationary scale, while its numerical algorithm was sub- 
sequently described in more detail in collaboration with 
Sastre p^. As in our spectral approach, the nonlinear 
evolution of the scalar fields is used to compute the source 
for the hij. However, |17j evolve the perturbation in po- 
sition space using a staggered leapfrog scheme, and we 
will refer to this algorithm as the leapfrog method. The 
implementation described in [171 imposes the transverse- 
traceless constraints at discrete times during the simula- 
tion, rather than at the source level as we do. 

The final algorithm is that of Dufaux, Bergman, 
Felder, Kofman and Uzan [TS], who develop a Green's 
function for the tensor portion of hij, which are evolved 
in Fourier space, as in our spectral method. We refer to 
this as the "Green's function" algorithm, and it is con- 
structed in an expanding background. The one explicit 
approximation in this algorithm is that it assumes the 
modes which are well inside the current Hubble horizon, 
or A: ^ aH, where k is the comoving wavenumber, a is 
the scale factor, and H is the Hubble parameter'^ 

The one case which has been publicly examined using 
all four algorithms is resonance following Xcj)'^ infiation. 
The box method (as implemented in [T^]) yields spectra 
with a significantly higher amplitude than those found 
with the spectral or Green's function methods, but that 
the latter methods agree very closely with each other. 
There is a distinction between methods which directly 
incorporate the expansion of spacetime, and we will see 
that the universe does grow significantly as resonance 
continues in the model - and expansion tends to di- 
lute the gravitational wave background. Moreover, care 
has to be used when computing the source term for the 
gravitational waves, as we will be taking the differences of 



derivatives computed on a spatial lattice, which requires 
the subtraction of terms of similar magnitude from one 
another, a well known scenario for inducing numerical 
error. The Green's function and spectral methods yield 
very similar results for the same models, the algorithms 
were developed and coded independently, and are con- 
ceptually distinct. We thus have considerable confidence 
that the gravitational wave background from preheating 
is being accurately evaluated.** 

In addition to the present approach, the direct produc- 
tion of metric perturbations via parametric resonance has 
been considered [191 EHl EJ EH EH EH, but this mech- 
anism is different from that considered here, and is less 
likely to produce a detectable signal. Further, Felder and 
Kofman have considered the "fragmentation" of the in- 
fiaton following reheating, which is closely related to the 
source of the gravitational waves discussed here [5S]. 

The structure of this paper is as follows. We begin by 
reviewing preheating, and following this we describe the 
numerical computation of the transverse-traceless com- 
ponent of Tfj^^ for a scalar field dominated universe, and 
the subsequent computation of the gravitational wave 
background, illustrating our approach with Xc/)'^ models. 
We then consider models driven by a quadratic infiaton 
potentials, a low scale hybrid model, and a massless quar- 
tic model with negative coupling. The last model has a 
very different resonance structure from the other scenar- 
ios that have been examined in this context, and we again 
find a substantial gravitational wave spectrum. Finally, 
we summarize the numerical issues we have identified, 
and discuss future prospects for work in this area. 



II. PARAMETRIC RESONANCE AND 
PREHEATING 

We now review parametric resonance of scalar fields 
[261 EH ES ES EHl [31]. Consider a toy model comprised 
of one classical infiaton field cf) and one (possibly massive) 
scalar field x 

- + Id^xd^'x - V{c^, x), (3) 

where 

1^(0, x) = V{cf>) + ^mlx' + \g^cb\\ (4) 

In this simple picture (j) is the infiaton, whereas x is a 
generic "matter" field. It is not necessary to specify V{4>) 
at this point; however, the specific form of (/) — x coupling 
term is important. We have additionally assumed that 



^ In principle, one could construct a Green's function without 
making this approximation, and in most cases of interest the 
relevant modes do satisfy k 2> aH. 



^ The leapfrog code produces a spectrum for the Ac/)'' background 
that is qualitatively similar to that of the Green's function and 
spectral methods. 
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there are no other x self interactions, i.e. A^x"*. There- 
fore X obeys 

x + iHx- ^y''x + mlx + g''<P'x = Q- (5) 

Adopting the Fourier transform convention 

~f{k,t)= j d\f{x,t)e^^^^-^, (6) 

and using fk as a shorthand for f{k, t) we expand x ™ 
terms of its momentum modes, 



Xk + iHxk 



+ m'+ .g202 ) xfe - 



(7) 



We can reduce this to a known analytic form by a) ignor- 
ing the expansion of the universe, so iJ = and k/a ^ k 
is simply the physical momentum, and b) setting rriy- — 0. 
Typically, the end of inflation is marked by oscillation of 
the inflation fleld about the minimum of its potential, 
and we will assume that this is described by a single 
quadratic term 



v{4>) 



1 



2 i2 



(8) 



With this potential i 
motion, so that 



' undergoes damped simple harmonic 



10 




FIG. 1: The Mathieu Stability/Instability chart. The imag- 
inary part of the Mathieu critical exponent is plotted, with 
darker colors corresponding to a larger imaginary component. 
Outside the heavy black lines the exponent is real-valued, and 
the corresponding solutions are strictly oscillatory. The diag- 
onal line corresponds to Aq — 2q. 



(f>{t) — <i>sin(TO^i), 



(9) 



where $ is a time-dependent amplitude, which varies 
slowly over a single cycle. Once we ignore the expan- 
sion of the universe, the friction term drops out of the 
equations of motion and is strictly constant. Substi- 
tuting into the mode equations and changing variables 
to 



fc2 



4771^ ' 2m? 
turns ([t]) into a Mathieu Equation [32]) 

Xfe + (^fc-29cos(2z))xfe = 



+ 2q, z = mt, (10) 



(11) 



where primes denote differentiation with respect to z. 
The solutions to a Mathieu equation are either stable 
(oscillatory) or unstable (exponential), depending on the 
relative values of q and Aq. Schematically, every solution 
to Mathieu's equation has two parts. 



Xk oc f{z)e 



:tif-LZ 



(12) 



where f{z) is periodic. When the Mathieu characteristic 
exponent /x has an imaginary part the solution has an 
exponentially growing mode. Figure [l] shows the values 
of S(/i) as a function of Ag and q. We see that ([To| 
requires Aq > 2g, so we are interested in the parameter 
values that lie to the left of the diagonal line in Figure [T| 



This ensures that the we find discrete resonance bands, 
which grow weaker as k grows larger. 

This process it not unique to m'^(j)'^ models. For in- 
stance, consider 



y(0) = -A0^ 



(13) 



In this case, the oscillations after the end of inflation are 
not sinusoidal; rather, they obey elliptic cosines [2H] . 



(X cn [ X — Xq, 



1 



(14) 



where a; is a dimcnsionless conformal time variable de- 
fined by 



' 9\ 1/4 

6 Am,, \ 



Vt. 



(15) 



In terms of these new variables, the mode equations for 
X are gH] 



Xk ■ 



where 0o is the initial amplitude of the (p field. This is 
a Lame equation |32j which has both exponential and 
oscillatory solutions, and can be related to the Mathieu 
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equation by writing 



V2 



27r\/2 
T 



n— 1 



-7r(n-l/2) 



_|_ g-7r(n-l/2) 



COS 



27r(2n- l)a; 



(17) 

In this expression, T is the period of the oscillation. The 
n — \ term has an amplitude of 0.9550 and dominates 
this series (see (28j for full analysis). If one neglects the 
remaining terms, the Mathieu equation is recovered. 

Generically, any periodic motion can be expanded in a 
sum of coherently oscillating terms. Each term in this ex- 
pansion will have a corresponding Mathieu equation that 
will lead to exponential excitation of particular momen- 
tum modes. However, the details of this expansion can 
vary significantly, and hence difi^erent models have very 
different resonant behavior. When we restore the expan- 
sion of the universe and allow mass terms to be non-zero 
along with any nonlinear interactions, the physics grows 
significantly more complicated as modes move in and out 
of resonance [26, 27, 28J. Further, this analysis neglects 
the evolution of (p; like the x field, the modes 4>k will be 
sourced by the g'^<\P"x^ coupling term. As these modes 
grow, (/) becomes inhomogeneous and the coherent oscil- 
lations will break down, due to backreaction from the 
created x particles. Ultimately, the only way to follow 
the full evolution is via numerical simulation. 

Standard treatments of parametric resonance describe 
the resonant pumping of momentum modes as particle 
creation. However, if we simply view and x as clas- 
sical Klein-Gordon fields, then the amplification of their 
Fourier modes makes the universe increasingly inhomo- 
geneous. The inhomogeneity in the fields leads to an 
inhomogeneous, time dependent energy density - which 
necessarily leads to the emission of gravitational radia- 
tion.^ 



III. GENERATION OF GRAVITATIONAL 
POWER 

A. The Power Spectrum 

The stress-energy tensor associated with gravitational 
radiation is given by |33j . 



T — 



1 



327rG 



\'<-ij,fj."- ,1// ) 



(18) 



and is specific to the transverse-traceless part of the met- 
ric perturbation. The brackets here denote a spatial av- 
erage and must be taken over a large enough volume in 
order to capture the contribution from long wavelength 
modes. In this work, we simply integrate over the full 



^ To be pedantic, this requires a non-zero quadrupole moment, but 
this is generically excited by an arbitrary perturbation. 



"grid" on which our numerical solutions are defined. The 
associated energy energy density is just the 00 compo- 
nent. 



Pgu 



1 



327rG 



(19) 



Recalling Parseval's theorem we can manipulate the spa- 
tial average to find 



1 

1? 



d3fc|F(fc)|2, 



(20) 



where L is the length of one side of the lattice. Thus 
the energy density can be expressed as an integral in 
momentum space. 



Pgw 



E 
E 
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2 

gig / dlnkk^\h,,,o{k)f 



(21) 
(22) 



or more usefully 
dp 



dink 



Lastly, this is 



dfl„ 



1 dp 



dink pcritdlnk 3H^L^ 



Ei^^^-.o(^)i' 



(23) 



(24) 



which we then insert into equation (20) of [12^^ to find, 

1/3 

.9 



o ^2 _ o u2d^gw{ae) ( go 



dink 



(25) 



where a,, is evaluated at the end of the simulation and 
go/ 9* is the ratio of number degrees of freedom today 
to the number degrees of freedom at matter/radiation 
equality and fi^ is the present day radiation energy den- 
sity. In this analysis we take go/g* = 1/100. 



B. Evolution of the Perturbation 



The goal our calculation will be to evaluate (19 1. By 



construction, the perturbation is free of scalar and vec- 
tor components, so it represents metric deformations due 
only to gravitational radiation. The perturbed Einstein's 
equations are 

G^,{t) + 5G^,{x., t) = S^G [f^,{t) + 5T^,{x, t)] (26) 



Please note a typographical error appears in I12| 
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where the background G^u and T^^^ obey Einstein's equa- 
tions for the unperturbed metric, 



and hence, 



(27) 



(28) 



One then obtains the equations of motion for the h 



%-2 



where 



2- I 

a , 



1. 



■%■] ^ ' I'll 
a a 



„2 '^ij ' 



(29) 



(30) 



and Slj is the transverse-traceless part of Sij . This can 
be extracted from Sij - as it can from any rank-2 tensor 
- by projecting onto the transverse plane and subtracting 
its trace [33 . Exphcitly, 



Sfj^ — PikSkiPij 



Pij (PlmSlm) (31) 



where the projection operator is defined by 



P — S 



(32) 



This procedure is implicitly required by the form of ( 29 1 



since one uses the transverse-traceless condition imposed 
on hij in order to derive this from a general metric per- 
turbation. Consequently, we must perform the same de- 
composition on our source in order to consistently evolve 
the metric perturbations. 

The hij obey (modified) wave equations, and in prin- 
ciple one could evolve D/i on the spatial lattice using 
the same numerical scheme employed to track the fields' 
evolution. In practice we found much better numerical 
stability when we wrote the hij in terms of their Fourier 
transforms, and solved the resulting ordinary differential 
equations for hij, 



G^,{k,t) = 8TrGf^,{k,t). 



(33) 



The fc = mode is our homogenous background, for 
which the corresponding component of hij necessarily 
vanishes. This is effectively a spectral algorithm for hij 
and we use a 4th order Runge-Kutta integrator to evolve 
the modes. At each point in space, we calculate Sij{x) 
explicitly from 



(34) 



and the full T^^ is obtained by summing the above ex- 
pression over all scalar fields. The nonvanishing compo- 



nents of Sij are 



di4'kdj4'k - 7;Sij [9m(?!)fe(9"Vfe] 



(35) 
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FIG. 2: From top to bottom: The scale factor, the variances 
of (j> (solid) and x (dotted), Sii" for a mode corresponding 
to |fc| ~ 1.4 X 10* Hz today, and the maximum height of the 
gravitational wave spectrum for this mode as a function of 
time for A<^'' inflation with A = 10"^* and g'/\ = 120. 



7 



We Fourier-transform this and construct 5/- (k), the 



source term for ( 29 1 . 

The field evolution was computed using LAT- 
TICEEASY [15 , which uses a staggered leapfrog inte- 
grator. The fields obey the Klein-Gordon equation in an 
expanding background, 



3°0, - 4v'<i 
a 



dV{(t)) 



0, 



(36) 



where the subscript i labels the fields. During resonance 
the gradient terms play a vital role, but in their absence 
we recover the familiar inflationary equation of motion. 
As usual, a is obtained from the Friedmann equations. 



a 



aV] =0. 



(37) 



Consequently, these simulations are performed in an ex- 
panding, rigid spacetime and we have ignored any backre- 
action from metric perturbations onto the field evolution. 
Numerical values for a(t) are normalized to unity at the 
beginning of our simulations, which start at the end of 
inflation. 

In Figure [2] we summarize the evolution of the grav- 
itational wave background for a V{(j)) = A0^/4 model, 
which is the de facto testbed for this type of calcu- 
lation, even though the underlying inflationary model 
is not consistent with recent CMB data. The specific 



model we solve is given by equation ( 13 1, with A = 10^^* 
and g^/X = 120. We explicitly check that the numeri- 
cal evolution of hij respects the transverse-traceless con- 
dition. Figures [3] and |4] explicity show representative 
Fourier components of the metric perturbation, along 
with explicit plots showing that that the numerical evolu- 
tion preserves the transversivity and tracelessness of hij . 
Note that the transversivity condition ^ takes the form 
kihii{k) + k2hi2{k) + k^hiz{k) in momentum-space. 

Unlike [TIT and fT^, we can compute the gravitational 
wave spectrum at any point during the simulation - a 
feature shared by the leapfrog [16] and Green's func- 
tion algorithms TH]. Figure [5] shows the evolution of 
the gravitational wave spectrum for V{4>) = X(f)'^ model 
with A = 10-^^ and g^/X = 120. The final amphtude of 
this spectrum is noticeably smaller than that computed 
using the "box" algorithm [T^]. However, this algorithm 
did not fully incorporate the expansion of space in its 
computation of the gravitational wave spectrum, and the 
universe undergoes significant expansion during the res- 
onant phase, and this expansion will have the effect of 
diluting the gravitational waves, and reducing their fi- 
nal amplitude. One can clearly see the "pumping" of 
the low frequency modes as resonance begins, and then 
the growth of the higher frequency modes as it continues. 
The evolution of hij is driven only by the inhomogeneities 
in the fields, and is not directly sourced by the poten- 
tial. The overall form of the field's power spectra, shown 
in Figure [6] closely resembles that of the gravitational 
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FIG. 3: The upper plot shows the magnitude, in arbitrary 
units, versus rescaled time of the diagonal components, /in 
(solid), hi2 (dotted) and his (dashed) for a mode correspond- 
ing to |fc| = 1.4 X 10* today for a Xij)'^ model. The lower plot 
shows k\h\i(k) -\- fc2/ii2(fc) -f kshiz{k), which demonstrates 
that perturbation is transverse via 



Our procedure applies to any source of stochastic grav- 
itational radiation where the signal is generated by large 
scale inhomogeneities, such as cosmological phase tran- 
sitions [3T and more specifically large scale bubble col- 
lisions [35j. Thus provided these processes can be sim- 
ulated numerically, this algorithm could be adapted to 
compute the gravitational radiation they generate. For 
the latter case, the crucial limitation is the resolution of 
the lattice itself: to accurately model the process of bub- 
ble nucleation till bubble collision, the distance between 
two neighbouring lattice grid points must be smaller than 
the nucleation size of the bubble, while the total lattice 
size must be larger than the bubble radius at percolation. 



C. Experimental Prospects 

Most theoretical discussions of gravitational wave spec- 
tra are expressed in terms of the ratio between the spec- 
tral energy density of the gravitational wave and the total 
energy density, as measured at the present day: 
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FIG. 4: The upper plot shows the magnitude, in arbitrary 
units, versus rescaled time of the diagonal components, hn 
(solid), h22 (dotted) and hss (dashed) for a mode correspond- 
ing to \k\ = 1.4 X 10* today for a X(fi'^ model. The lower plot 
shows their sum. 
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FIG. 6: The evolution of the power spectra of the fields for 
both (j) (above) and x (below) . This is done for the case where 
A = 10"^" and g^/X = 120 on a 128^ lattice. 
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FIG. 5: The evolution of the gravitational wave spectrum for 
a quartic inflation model, where A = 10^^'' and g^/\ = 120. 



Experimental bounds are commonly expressed in terms 
of the spectral density, Sh{f), which is an ensemble av- 
erage over the Fourier amplitudes of the plane-wave so- 
lutions /i^(/, f2) averaged over all directions [37]. That 
is, a transverse-traceless metric perturbation can be de- 
composed into plane wave solutions, each carrying two 



polarizations, in all directions, i.e 

''^1= E r df f dnhA{f,n)e-^^'f'efj{n) (39) 

Here, efj is a tensor that identifies the polarizations rel- 
ative to the coordinates and and f2 is a directional 
unit vector with differential, dfl = d cos 6d(j). From this, 

s{f - f)SAA'\sh{f) ^ J dndn' (h\{f,n)hA'{f,n')y 

(40) 

More usefully, this quantity can be directly related to 
flgw through 

47r2 

= ^f^l^if)- (41) 

We will typically find signals where ^lgw{f) is, at most, 
0(10""'^"'^). Consequently, in order to detect this signal 
we need a an observatory that can detect spectral den- 
sity on the order of Sk ^ lO-^^h'^Bz'^ at 100 Hz or 
Sh ^ lO-^'i/i^Hz"^ at 10-2 Hz. For comparison, the 
minimum spectral density detectable by Advaced LIGO 
is approximately 10"''^ Hz"^ at 100 Hz |3S], and LISA is 
10-*° Hz"^ at 10-2 Hz ISg. 
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IV. COSMOLOGICAL MODELS 



A. Quadratic Inflation 



30 



Having established our overall methodology, we now 
compute representative spectra for several explicit cos- 
mological models. In Section |IIIB| we presented results 
for a chaotic model with a quartic potential, so as to best 
compare our calculations to those of other groups. We 
now consider the potential 



20 



(42) 



100 



which was also treated in [12", 114]. The overall evolu- 
tion of this system is summarized in Figure |7] while the 
spectrum is similar to that obtained using the box algo- 
rithm [T?]. In this case the agreement between the box 
algorithm and the spectral code is considerably better, 
and we see that the box code actually underestimates 
the maximal power in gravitational waves. Note that the 
total growth of the universe during resonance is smaller 
in this case than it was during quartic inflation - and 
since the box algorithm effectively compute hij in flat 
space, the discrepancy introduced by this approximation 
is accordingly reduced. 



B. Low-Scale Inflation 

In inflationary models with a single free parameter in 
the potential, the scale of inflation is typically fixed by 
matching to the observed amplitude of the scalar pertur- 
bation spectrum. Consequently, for the quadratic poten- 
tial we discussed above, m is not a free parameter. Now 
consider a general hybrid inflation model [TTl 301 SI] , 



(Af2 - \a^) 



2\2 



2, ,2 



4A 2 ^ 

(43) 

So long as i/i is large, cr = is a stable minimum and 
the energy density is dominated by the potential energy 
associated with (p. When (j) « M/h, cr = is no longer a 
minimum and a quickly settles into one of the two new 
minima, a — a process known as the "water- 

fall" phase transition. We assume that the field a settles 
coherently into its vacuum expectation value and there 
is no a particle production. This is consistent with the 
assumption that m <^ M . Assuming, with no loss of gen- 
erality, that the field settles into its positive minimum, 
acquires the following effective potential energy. 



v{4>) - - + 



A 



or, since we have assumed a small m. 
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(44) 
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FIG. 7: From top to bottom: The scale factor, the variances 
of cj> (dotted) and x (solid), for a mode corresponding 

to |fc| = 5.4 X 10* Hz today, and the maximum height of the 
gravitational wave spectrum as a function of time. This is for 
the m = 10~^mp and g^mli/m^ = 2.5 x 10^ model. 
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In this regime, we recover a chaotic inflation model where 
4> acquires an effective mass, rn^ff = h?M^ /\. This sys- 
tem was studied in detail in [T^ by the present authors. 
In particular, we used this model to demonstrate that 
the spectrum of gravitational radiation produced during 
resonance has the scaling properties conjectured in [l2] . 
However, this model has the disadvantage of needing ex- 
ceptionally small parameter values in order to sure that 
the structure of resonance does not depend directly on 
the inflationary scale. 

C. General Hybrid Potentials 

Now consider the potential 

n0,x)=^0^+^X^+|0V, (46) 

where we allow g to aquire either sign. In principle |30j . 
one can include mass terms for both the (j) and the % 
flelds. We omit them here for simplicity, and we one 
must also ensure that if g < 0, 1^(0, x) ^ everywhere, 
or 

^ > 1. (47) 
9 

When (j) is large the x fi^ld has a substantial effective 
mass and remains at x = 0. If g > 0, the effective mass 
of X only vanishes at the origin of the potential. This 
suppresses the production of x particles and leads to the 
slow and inefficient production of (p particles. Conversely, 
for g < Q the x held picks up a tachyonic mass when 
(fp' is small, there are rapid oscillations in both and 
X, and particle production proceeds at a dramatic pace. 
While Figure [T| strictly applies to the Mathieu equation, 
this system effectively has a negative q and the resonant 
modes thus live in the lower portion of the plot, where 
the instability regions are broad, and the imaginary part 
of the critical exponent is large. Consequently, resonance 
occurs extremely promptly, and the resulting dynamics 
differ dramatically from those of all the other models 
considered here. In particular, the universe grows by a 
factor of ~ 3 as the vast majority of the gravitational 
wave power is generated. Further, this model produces 
(j) quanta, and the corresponding variance has a sharp 
peak, rather than the broad "hump" seen in the other 
cases we examine. 

We have chosen to work with the same parameters used 
in the simulations of [3D], = 10~^^, A,^ = 10"'', and 
\g\ = 10"^*^, which provides a cross-check on our evalu- 
ation of the fleld evolution. We summarize the resonant 
epochs in these models in Figures |8] and |9] for positive 
and negative g respectively. The variances in Figure [8] 
and Figure [9] are consistent with Figures 4 and 5 of [5U] . 
With a negative coupling constant the fields quickly be- 
come inhomogeneous. Conversely, with a positive cou- 
pling case the x field initially becomes inhomogeneous, 
with the (p following considerably later, leading to two 
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FIG. 8: From top to bottom: The scale factor (normalized 
to one at the beginning of the simluation), the variances of cj) 
(solid) and x (dotted), and one component, Si^ for a mode 
corresponding to |fc| « 5.4 x 10^ Hz today, and the maximum 
height of the gravitational wave spectrum as a function of 
time. This is for the g — 10~^° case. 
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FIG. 9: From top to bottom: The scale factor (normalized 
to one at the beginning of the simluation), the variances of 
(solid) and x (dotted), and one component, S^i" for a mode 
corresponding to |fc| « 5.2 x 10^ Hz today, and the maximum 
height of the gravitational wave spectrum as a function of 
time. This is for the g = —10"^" case. 
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FIG. 10: The overall spectrum (top panel) and the individual 
spectra (bottom panel) for the simulation shown in Figure [s] 
In the lower plot, the dotted line shows gravitational waves 
due to (/)-particle production while the solid line denotes the 
contribution from the x field. 
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FIG. 11: The spectrum of gravitational waves generated dur- 
ing preheating with a negative coupling, with the same pa- 
rameters as the simulation shown in Figure |9] 



superimposed sets of gravitational waves, as can be seen 
from the "crossing" of the two variances in Figure [8] 
We display the gravitational power spectra of these 

In the positive coupling 



models in Figures 10 and 11 



case, we separate the contribution from the two fields, 
and show that these sum together to give a spectrum 
with a pronounced peak. Conversely, when the coupling 
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is negative, particle production occurs very rapidly, and 
the associated gravitational wave spectrum again has a 
substantial amplitude, despite the very different reso- 
nance dynamics. We plan to examine this system more 
carefully in a future publication. 



NUMERICAL ISSUES 



The source term for hij contains differences of combi- 
nations of derivatives which arc evaluated on a discretized 
lattice, there is ample opportunity for numerical "noise" 
to contaminate these simulations. Consequently, we have 
to take care that the spatial extent of our lattice (hence- 
forth the "box") is small enough to ensure that highest 
Fourier modes that contribute to the gravitational wave 
spectrum are well resolved. 



In Figure 12 we show the consequences of making the 
box too large, and thus losing the fine resolution needed 
to properly resolve the peak. The specific mode l he re is 
V{(j)) = A0'*/4 with g'^/X = 120 as in Section ^ We 
see that the discretization noise typically adds spurious 
power to the spectrum. This issue does not appear to 
explain the discrepancy between the box algorithm and 
the spectral algorithm discussed in Section |III| - the lat- 
tice used in [12^ is large enough to resolve the peak, and 
the "box" algorithm yields a much better match to the 
results found for a quadratic potential. However, it is 
noteworthy that the universe grows significantly during 
the reheating phase following infiation with a quartic po- 
tential, and this growth would tend to dilute the gravita- 
tional waves as they are produced - and the box method 
does not properly account for the growth of the universe 
in its computational of the gravitational waves. Con- 
sequently, we conclude that it is safest to solve the hij 
self-consistently in an expanding background when evalu- 
ating the gravitational wave signal generated during pre- 
heating. Further the computational cost of doing so is 
relatively small. 

In Figure[l3]we compare results from our spectral code 
to those obtained using the Green's function algorithm 
of Dufaux et al. |18j . and find that the two methods 
overlap exceptionally well. For the "standard" 128"^ lat- 
tice used in the calculations in this paper, the agreement 
between the two methods is good but not perfect. We 
also present data from 256'^ and 512^ runs which agree 
very closely with one another, and are presumably close 
to the continuum limit. ^ Moreover, these results also 
overlap closely with the Green's function result, which 
provides a crucial cross-check oft the algorithms used by 
ourselves and the authors of J^. Note that both codes 
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^ Hal Finkel collaborated with us on writing the code used for the 
256^ and 512^ runs, and we thank Gary Felder with providing 
us with the raw data used to plot the Green's function spectrum 
in Figure [Tsl 



FIG. 12: This shows a test of the integration step size. The 
dotted lines show simulations, with varying box sizes, where 
the metric perturbations are evolved at every timestep, the 
solid line shows the case where the timestep for the metric 
perturbation is increased. 



use the same transfer function to shift these spectra into 
their presently observable values. Moreover, the results 
of |T7] appear to be in good qualitative agreement with 
those presented here. 

The numerical simulations here are conceptually 
straightforward - we are, after all, simply evolving a set 
of second order partial differential equations, but expe- 
rience has shown that their actual implementation can 
require some finesse. We now have solid agreement be- 
tween two independent codes - the discrepancies between 
the two algorithms shown in Figure[T3]are of the same or- 
der as the changes induced in the spectrum by modifying 
the size of the underlying lattice. One can thus be confi- 
dent that the gravitational wave signal generated during 
preheating is being accurately evaluated, and these re- 
sults can be safely used to assess observational strategies 
for detecting a stochastic background. 

We do note, however, that almost all the different nu- 
merical approaches summarized in the introduction have 
used a staggered leapfrog algorithm to evolve the fields, 
and it would be worth exploring the consequences of us- 
ing different algorithms for this part of the calculation, 
particularly those with better than second order accuracy 
in their spatial differencing. In particular, since we are 
using a spectral algorithm for the metric perturbations, 
it is a natural extension of our approach to implement 
a spectral solver for the field evolution as well. This is 
somewhat complicated since the couplings between the 
fields render these equations nonlinear, but the applica- 
tion of spectral methods to nonlinear systems is a well 
understood problem. 

As presently implemented, our code uses a good deal 
of memory ~ each field requires two "copies" of the grid 
(the field, and its first time derivative), plus 12 more 
to evolve the Fourier modes of the 6 hij , and a further 
handful for computing source terms and Fourier trans- 
forms. If we assume isotropy (which is a very reasonable 
expectation) , equation 24 is a sum over three equivalent 
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FIG. 13: We show the resuhs of simulation a quartic inflation with jX — 120 using both our spectral code, and the Green's 
function algorithm. This simulation matches Figure 4 of 18 , and the Green's function results are shown with a dotted line. 
We show results from the spectral code for three different lattice sizes - 128^ (solid), 256^ (dashed) and 512^ (dot-dashed). 
The last two lines overlap, so we are presumably approaching the continuum limit, and we also have excellent agreement with 
the results of [18] , in which the gravitational wave spectrum was obtained using a very different numerical algorithm. 



diagonal terms and three equivalent off-diagonal terms. 
Consequently, we can reduce the labor and storage in- 
volved in computing the hij by a factor of three if we 
assume that h\\ — /i22 = /133 and h\2 = /ii3 = /i23, and 
this works as expected in practice. We could further im- 
prove matters by only evolving a subset of the Fourier 
modes of the hij (a strategy employed by [18]) - the only 
downside to these approaches is that, if carried too far, 
the spectrum acquires significant stochastic noise, since 
it is being computed from a small sample of total set of 
modes. 

The evolution code of the fields will have some charac- 
teristic timestep dt^ whose optimum value is a function 
of the specific dynamical system, and the lattice spac- 
ing dx. Conversely, most of the modes of the hij can be 
evolved with a significantly larger timestep, so one has 
the option of updating the hij evolution on a different 
timescale to that of the fields, further improving perfor- 
mance.^ In this analysis, which is essentially a "proof of 



concept" , we have not aggressively pursued these opti- 
mizations but we do plan to employ them in our future 
work. In our current implementation, the evaluation of 
hij takes approximately 90% of the overall runtime, so 
any optimization to this calculation will have substantial 
benefits. 



VI. DISCUSSION 

This paper describes the spectral algorithm for comput- 
ing a stochastic background of gravitational waves gen- 
erated during an epoch in which the universe is highly 
inhomogeneous on small scales. We have focussed on 
the spectrum generated at the end of infiation, during an 
era of preheating and/or parametric resonance. However, 
this algorithm is applicable to any epoch in which inho- 
mogeneities source gravitational radiation, including first 
order phase transitions (34j , bubble collisions (35j i36j , or 
MHD turbulence [i^H5]. 



This strategy is particularly suited to cases where one evolves the 
Fourier modes of the hij since the mode equations are indepen- 
dent at linear order in perturbation theory. If one chooses to solve 
Ohij on a spatial grid the timestep will need to be commensurate 



with the lattice spacing in order to satisfy the Courant condition, 
unless one "downsamples" the grid for the metric terms, relative 
to that of the fields. 
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We began this paper with a review of other recent ap- 
proaches to this problem, highhghting their relative mer- 
its. In addition to the previous work of the current au- 
thors [T21[T3], two other independent codes have been re- 
cently been developed to address this problem [THl [T71 [TS] 
and we explicitly compare the results of our spectral algo- 
rithm that the Green's function algorithm of TH], show- 
ing that the difference between the computed spectra is of 
the same order as that introduced by the finite resolution 
of the underlying spatial lattice. This is a welcome devel- 
opment, as it validates both algorithms and their explicit 
implementations. Conversely, algorithms for evaluating 
the gravitational wave spectrum which do not account 
for the expansion of the universe do not necessarily over- 
lap with the results obtained with our spectral algorithm, 
and should not be used if the universe expands by a large 
(0(10) or more) factor during the resonant phase. 

The overall scaling properties of any stochastic spec- 
trum produced as inflation ends and the universe reheats 
were first conjectured in [T5]. These claims were con- 
firmed explicitly for a simple resonant model in [14j , using 
the algorithm we describe here. In this paper we extend 
our analysis to a class of hybrid models and find spectra 
that are qualitatively different from those seen previous 
in this context. In particular, models with an "inverted" 
coupling, first discussed in [30 , have a particularly strong 
period of resonance. This is a significant development for 
several reasons - firstly, resonance here is substantially 
different from that of the chaotic models investigated pre- 
viously. Secondly, this model provides a more realistic ex- 
ample of resonance with an arbitrary inflationary scale. 
We plan to investigate its properties more carefully in a 
future paper, in order to better understand the correla- 



tion between they dynamics of resonance and preheating 
and the resulting gravitational wave background. 

In this paper, for the sake of computational conve- 
nience most of our calculations were performed on 128^ 
grids, which easily fits on a single computational node. 
However, our spectral code has been extended to a clus- 
ter environment, and which allows simulations on much 
larger grids. We plan to use this code to investigate the 
full dynamic range of the spectrum, especially in mod- 
els where the post-inflationary Hubble scale is not vastly 
larger than the scales which undergo resonance. Finally, 
we also plan to apply our algorithm to other situations 
where signiflcant gravitational wave backgrounds can be 
generated by small-scale inhomogeneities in the primor- 
dial universe. 
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